This tutorial is part of viskex.

import typing

import firedrake
import numpy as np
import ufl

/home/firedrake/firedrake/lib/python3.10/site-packages/pytools/__init__.py:2449: UserWarning: unable to find git revision
  warn("unable to find git revision")
firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance

import viskex

Generate and plot meshes

interval = firedrake.UnitIntervalMesh(6)
square_tria = firedrake.UnitSquareMesh(3, 3, quadrilateral=False, diagonal="right")
square_quad = firedrake.UnitSquareMesh(3, 3, quadrilateral=True, diagonal="right")
cube_tetra = firedrake.UnitCubeMesh(3, 3, 3, hexahedral=False)
cube_hexa = firedrake.UnitCubeMesh(3, 3, 3, hexahedral=True)

for mesh in (interval, square_tria, square_quad, cube_tetra, cube_hexa):
    mesh.topology_dm.removeLabel(firedrake.cython.dmcommon.FACE_SETS_LABEL)
del mesh

viskex.firedrake.plot_mesh(interval)

viskex.firedrake.plot_mesh(interval, dim=1)

viskex.firedrake.plot_mesh(interval, dim=0)

viskex.firedrake.plot_mesh(square_tria)

viskex.firedrake.plot_mesh(square_tria, dim=2)

viskex.firedrake.plot_mesh(square_tria, dim=1)

viskex.firedrake.plot_mesh(square_tria, dim=0)

viskex.firedrake.plot_mesh(square_quad)

viskex.firedrake.plot_mesh(square_quad, dim=2)

viskex.firedrake.plot_mesh(square_quad, dim=1)

viskex.firedrake.plot_mesh(square_quad, dim=0)

viskex.firedrake.plot_mesh(cube_tetra)

viskex.firedrake.plot_mesh(cube_tetra, dim=3)

viskex.firedrake.plot_mesh(cube_tetra, dim=2)

viskex.firedrake.plot_mesh(cube_tetra, dim=1)

viskex.firedrake.plot_mesh(cube_tetra, dim=0)

viskex.firedrake.plot_mesh(cube_hexa)

viskex.firedrake.plot_mesh(cube_hexa, dim=3)

viskex.firedrake.plot_mesh(cube_hexa, dim=2)

viskex.firedrake.plot_mesh(cube_hexa, dim=1)

viskex.firedrake.plot_mesh(cube_hexa, dim=0)

Generate and plot subdomains

def mark_subdomains(mesh: firedrake.MeshGeometry) -> firedrake.MeshGeometry:  # type: ignore[no-any-unimported]
    """Mark left and right subdomains in a given mesh with values 1 and 2, respectively."""
    cellname = mesh.ufl_cell().cellname()
    if cellname in ("interval", "triangle", "tetrahedron"):
        subdomains_function_space = firedrake.FunctionSpace(mesh, "DP", 0)
    elif cellname in ("quadrilateral", "hexahedron"):
        subdomains_function_space = firedrake.FunctionSpace(mesh, "DQ", 0)
    else:
        raise RuntimeError("Invalid cellname")
    x = firedrake.SpatialCoordinate(mesh)
    left_subdomain = firedrake.Function(subdomains_function_space).interpolate(
        firedrake.conditional(x[0] <= 1.0 / 3.0, 1.0, 0.0))
    right_subdomain = firedrake.Function(subdomains_function_space).interpolate(
        firedrake.conditional(x[0] >= 2.0 / 3.0, 1.0, 0.0))
    mesh_with_subdomains = firedrake.RelabeledMesh(mesh, [left_subdomain, right_subdomain], [1, 2])
    mesh_with_subdomains.init()
    return mesh_with_subdomains

interval_with_subdomains = mark_subdomains(interval)
square_tria_with_subdomains = mark_subdomains(square_tria)
square_quad_with_subdomains = mark_subdomains(square_quad)
cube_tetra_with_subdomains = mark_subdomains(cube_tetra)
cube_hexa_with_subdomains = mark_subdomains(cube_hexa)

viskex.firedrake.plot_mesh_entities(
    interval_with_subdomains, 1, "subdomains", interval_with_subdomains.cell_subset(2).indices
)

viskex.firedrake.plot_mesh_entities(
    interval_with_subdomains, 1, "subdomains", interval_with_subdomains.cell_subset(2).indices,
    2 * np.ones_like(interval_with_subdomains.cell_subset(2).indices)
)

viskex.firedrake.plot_mesh_sets(interval_with_subdomains, 1, "subdomains")

viskex.firedrake.plot_mesh_entities(
    square_tria_with_subdomains, 2, "subdomains", square_tria_with_subdomains.cell_subset(2).indices
)

viskex.firedrake.plot_mesh_entities(
    square_tria_with_subdomains, 2, "subdomains", square_tria_with_subdomains.cell_subset(2).indices,
    2 * np.ones_like(square_tria_with_subdomains.cell_subset(2).indices)
)

viskex.firedrake.plot_mesh_sets(square_tria_with_subdomains, 2, "subdomains")

viskex.firedrake.plot_mesh_sets(square_quad_with_subdomains, 2, "subdomains")

viskex.firedrake.plot_mesh_sets(cube_tetra_with_subdomains, 3, "subdomains")

viskex.firedrake.plot_mesh_sets(cube_hexa_with_subdomains, 3, "subdomains")

Generate and plot boundaries

def mark_boundaries(mesh: firedrake.MeshGeometry) -> firedrake.MeshGeometry:  # type: ignore[no-any-unimported]
    """
    Mark internal and boundary facets in a given mesh with four different values.

    Internal facets of left and right subdomains are associated with values 1 and 2, respectively.
    Furthermore, boundary facets on the left and right boundaries are associated with values 3 and 4,
    respectively.
    """
    cellname = mesh.ufl_cell().cellname()
    if cellname in ("interval", ):
        boundaries_function_space = firedrake.FunctionSpace(mesh, "P", 1)
    elif cellname in ("triangle", "quadrilateral", "tetrahedron", "hexahedron"):
        boundaries_function_space = firedrake.FunctionSpace(mesh, "HDiv Trace", 0)
    else:
        raise RuntimeError("Invalid cellname")
    x = firedrake.SpatialCoordinate(mesh)
    left_boundary = firedrake.Function(boundaries_function_space).interpolate(
        firedrake.conditional(abs(x[0] - 0.) < np.sqrt(np.finfo(float).eps), 1.0, 0.0))
    right_boundary = firedrake.Function(boundaries_function_space).interpolate(
        firedrake.conditional(abs(x[0] - 1.) < np.sqrt(np.finfo(float).eps), 1.0, 0.0))
    left_subdomain = firedrake.Function(boundaries_function_space).interpolate(
        firedrake.conditional(
            firedrake.And(x[0] <= 1.0 / 3.0, abs(x[0] - 0.) > np.sqrt(np.finfo(float).eps)), 1.0, 0.0))
    right_subdomain = firedrake.Function(boundaries_function_space).interpolate(
        firedrake.conditional(
            firedrake.And(x[0] >= 2.0 / 3.0, abs(x[0] - 1.) > np.sqrt(np.finfo(float).eps)), 1.0, 0.0))
    mesh_with_boundaries = firedrake.RelabeledMesh(
        mesh, [left_boundary, right_boundary, left_subdomain, right_subdomain], [3, 4, 1, 2])
    mesh_with_boundaries.init()
    return mesh_with_boundaries

interval_with_boundaries = mark_boundaries(interval)
square_tria_with_boundaries = mark_boundaries(square_tria)
square_quad_with_boundaries = mark_boundaries(square_quad)
cube_tetra_with_boundaries = mark_boundaries(cube_tetra)
# cube_hexa_with_boundaries = mark_boundaries(cube_hexa)  # HDiv not implemented

viskex.firedrake.plot_mesh_entities(
    interval_with_boundaries, 0, "boundaries",
    interval_with_boundaries.exterior_facets.measure_set("exterior_facet", 4).indices
)

viskex.firedrake.plot_mesh_entities(
    interval_with_boundaries, 0, "boundaries",
    interval_with_boundaries.exterior_facets.measure_set("exterior_facet", "everywhere").size
    + interval_with_boundaries.interior_facets.measure_set("interior_facet", 2).indices
)

viskex.firedrake.plot_mesh_sets(interval_with_boundaries, 0, "boundaries")

viskex.firedrake.plot_mesh_entities(
    square_tria_with_boundaries, 1, "boundaries",
    square_tria_with_boundaries.exterior_facets.measure_set("exterior_facet", 4).indices
)

viskex.firedrake.plot_mesh_entities(
    square_tria_with_boundaries, 1, "boundaries",
    square_tria_with_boundaries.exterior_facets.measure_set("exterior_facet", "everywhere").size
    + square_tria_with_boundaries.interior_facets.measure_set("interior_facet", 2).indices
)

viskex.firedrake.plot_mesh_sets(square_tria_with_boundaries, 1, "boundaries")

viskex.firedrake.plot_mesh_sets(square_quad_with_boundaries, 1, "boundaries")

viskex.firedrake.plot_mesh_sets(cube_tetra_with_boundaries, 2, "boundaries")

Interpolate and plot scalar functions

def prepare_scalar_field_cases(  # type: ignore[no-any-unimported]
    mesh: firedrake.Mesh,
    expression: typing.Callable[[ufl.core.expr.Expr], ufl.core.expr.Expr]
) -> typing.Tuple[firedrake.Function, typing.Tuple[ufl.core.expr.Expr, ufl.FunctionSpace]]:
    """Prepare scalar field cases."""
    scalar_function_space = firedrake.FunctionSpace(mesh, "CG", 2)
    scalar_field_ufl = expression(ufl.SpatialCoordinate(mesh))
    scalar_field = firedrake.interpolate(scalar_field_ufl, scalar_function_space)
    return scalar_field, (scalar_field_ufl, scalar_function_space)

interval_scalar_field, interval_scalar_field_ufl = prepare_scalar_field_cases(
    interval, lambda x: x[0]**3)
square_tria_scalar_field, square_tria_scalar_field_ufl = prepare_scalar_field_cases(
    square_tria, lambda x: x[0]**3 + x[1]**2)
square_quad_scalar_field, square_quad_scalar_field_ufl = prepare_scalar_field_cases(
    square_quad, lambda x: x[0]**3 + x[1]**2)
cube_tetra_scalar_field, cube_tetra_scalar_field_ufl = prepare_scalar_field_cases(
    cube_tetra, lambda x: x[0]**3 + x[1]**2 + x[2]**4)
cube_hexa_scalar_field, cube_hexa_scalar_field_ufl = prepare_scalar_field_cases(
    cube_hexa, lambda x: x[0]**3 + x[1]**2 + x[2]**4)

viskex.firedrake.plot_scalar_field(interval_scalar_field, "scalar")

viskex.firedrake.plot_scalar_field(interval_scalar_field_ufl, "scalar")

viskex.firedrake.plot_scalar_field(square_tria_scalar_field, "scalar")

viskex.firedrake.plot_scalar_field(square_tria_scalar_field, "scalar", warp_factor=1.0)

viskex.firedrake.plot_scalar_field(square_tria_scalar_field_ufl, "scalar")

viskex.firedrake.plot_scalar_field(square_quad_scalar_field, "scalar")

viskex.firedrake.plot_scalar_field(square_quad_scalar_field, "scalar", warp_factor=1.0)

viskex.firedrake.plot_scalar_field(square_quad_scalar_field_ufl, "scalar")

viskex.firedrake.plot_scalar_field(cube_tetra_scalar_field, "scalar")

viskex.firedrake.plot_scalar_field(cube_tetra_scalar_field_ufl, "scalar")

viskex.firedrake.plot_scalar_field(cube_hexa_scalar_field, "scalar")

viskex.firedrake.plot_scalar_field(cube_hexa_scalar_field_ufl, "scalar")

Interpolate and plot vector functions

def prepare_vector_field_cases(  # type: ignore[no-any-unimported]
    mesh: firedrake.Mesh,
    expression: typing.Callable[[ufl.core.expr.Expr], ufl.core.expr.Expr]
) -> typing.Tuple[firedrake.Function, typing.Tuple[ufl.core.expr.Expr, ufl.FunctionSpace]]:
    """Prepare vector field cases."""
    vector_function_space = firedrake.VectorFunctionSpace(mesh, "CG", 2)
    vector_field_ufl = ufl.as_vector(expression(ufl.SpatialCoordinate(mesh)))
    vector_field = firedrake.interpolate(vector_field_ufl, vector_function_space)
    return vector_field, (vector_field_ufl, vector_function_space)

square_tria_vector_field, square_tria_vector_field_ufl = prepare_vector_field_cases(
    square_tria, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
square_quad_vector_field, square_quad_vector_field_ufl = prepare_vector_field_cases(
    square_quad, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
cube_tetra_vector_field, cube_tetra_vector_field_ufl = prepare_vector_field_cases(
    cube_tetra, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))
cube_hexa_vector_field, cube_hexa_vector_field_ufl = prepare_vector_field_cases(
    cube_hexa, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))

viskex.firedrake.plot_vector_field(square_tria_vector_field, "vector")

viskex.firedrake.plot_vector_field(square_tria_vector_field, "vector", glyph_factor=0.1)

viskex.firedrake.plot_vector_field(square_tria_vector_field, "vector", warp_factor=1.0)

viskex.firedrake.plot_vector_field(square_tria_vector_field_ufl, "vector")

viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector")

viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector", glyph_factor=0.1)

viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector", glyph_factor=0.1)

viskex.firedrake.plot_vector_field(square_quad_vector_field, "vector", warp_factor=1.0)

viskex.firedrake.plot_vector_field(cube_tetra_vector_field, "vector")

viskex.firedrake.plot_vector_field(cube_tetra_vector_field, "vector", glyph_factor=0.1)

viskex.firedrake.plot_vector_field(cube_tetra_vector_field, "vector", warp_factor=1.0)

viskex.firedrake.plot_vector_field(cube_hexa_vector_field, "vector")

viskex.firedrake.plot_vector_field(cube_hexa_vector_field, "vector", glyph_factor=0.1)

viskex.firedrake.plot_vector_field(cube_hexa_vector_field, "vector", warp_factor=1.0)